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A couple of orthogonal coordinates for dipole geometry are proposed for numerical simulations 
of plasma geophysics in the Earth's dipole magnetic field. These coordinates have proper metric 
profiles along field lines in contrast to the standard dipole coordinate system that is commonly used 
in analytical studies for dipole geometry. 



I. INTRODUCTION 

In the study of plasma geophysics, an orthogonal coordinate system defined by a dipole field is commonly used 
because of the Earth's dipole magnetic field B^. The standard dipole coordinate system (/i, x, 0) is defined through 
the spherical coordinates (r, 6, cj)) as 
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where r is length from Earth's center, normalized by its radius li?e, 9 is colatitude, and <j} is the longitude. The 
coordinate /i is a potential function of a dipole field, oc V/i, and constant-x curves in a meridian plane, (j) — const., 
denote dipole filed lines. 

Since (/i, x, <j)) is an orthogonal system, their metric terms are simply given by 
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with 



6(61) = + 3 cos2 9, (5) 



and the length clement ds is given by ds^ — ds^ + ds^ + ds^ with ds^ — d^,, ds^ = dx, dSc/> = ft-^ dcj). Given 
these metric terms, it is straightforward to discretize any differential operator such as the divergence of a vector v 
that is denoted by components {u^, w^, v^} in the dipole coordinates as 
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by, for example, a finite difference method in the computational (/i, x, (f) space. 

The above standard dipole coordinates is convenient and certainly appropriate for analytical studies in which the 
Earth's dipolar field plays central roles. It also works as a base coordinates for the node and cell generation of the 
finite element method in the dipole geometry However, when one tries to use other numerical methods in which 

analytical expression of the metric terms are important for preserving numerical simplicity and accuracy, as in the 
case of the finite difference method, the standard dipole coordinate x, 0) cannot be used in its original form since 
the metric changes intensely along the field lines. 

It should be noted that oc jB^j^^ from the above definitions, which means that /i^ is roughly proportional to 
r^. Therefore, the metric /i^ at r = 1 is O(IO^) smaller than that at r = 10. Fig.^a) shows the profile along a 
field line starting from 70°N as a function of fx. (We suppose that the north pole is located in 6' = in this note.) 
This field fine goes through the equator (/i = 0) at r = 8.55. Note the sharp peak in Fig.^a) at the equator. 

When one uses the finite central difference method, the grid spacing along the field line is given by As^ = Afi. 
Fig.|21shows grid point distribution in the standard dipole coordinates. The grid size in the figure is x — 101 x 10. 
(101 grids along each field line and 10 grids in the perpendicular direction.) The starting points of the field lines are 
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between 65°N and 70°N at r = 1. All the grid points are shown in the figure without any skip. It is clearly seen that 
the resolution near the equator is so poor that any numerical simulation on this grid system is impractical. For the 
field line starting from 70°N at r = 1, the metric ft.^ on the equator is about 1160 times larger than that on r = 1. 
Also note that the imbalance of the grid spacings between the near Earth and the near equatorial regions along the 
field lines causes unnecessarily severe restriction on the Courant-Friedrichs-Lewy condition in explicit time integration 
schemes. 

For some numerical simulations such as the magnetosphere-ionosphere coupling, it is certainly desirable to use a 
grid system that has a natural grid convergence near the ionosphere to resolve fine structures near the coupling region, 
but the three-orders of magnitude is obviously too much. This is especially serious when one tries to simulate some 
phenomena in which relatively high resolution near the equator is required. An example of such simulation is the 
auroral arc formation by the feedback instability driven by vortex flow in the equator |llL Il2j | . 

A trivial way to avoid the poor resolution problem of the standard dipole coordinates near the equator is to place 
the computational grid points along the /x space in a nonuniform way. In this case, the metric factors have to be 
numerically calculated. For example, Lee and Lysak determined the grid spacing due to the local Alfven wave 
speed. The same approach was adopted in Budnik et al. [l|. However, this method injures the generality of the dipole 
coordinates as well as its analytical nature. 

If one prefers to fully numerical methods, refer to Proehl et al. in which a general algorithm to construct grid 
points along an arbitrarily given magnetic field, including the dipole, is presented. In contrast to that approach, 
we propose in this note analytical as well as simple coordinate transformations of fi that lead to practical metric 
distributions along the field line. 

II. TRANSFORMATION FORMULA OF THE DIPOLE COORDINATES 

Before we go into the description on the modified dipole coordinates defined by the coordinate transformation of the 
standard dipole coordinates, we derive analytical expressions of the inverse transformation from the standard dipole 
coordinates {/x, x, (p) into the spherical coordinates (r, 9, (j)) since we could not find these expressions in the literature 
and they can be directly applied to the modified dipole coordinates described later. 

Eliminating r from in eq. with subsidiary variables u and C defined as 

u^sinH, C=Wx')', (7) 

we get a fourth order equation of u: 

(^" + ^-1 = 0. (8) 



The unique solution of eq. Q for positive real u is 




(9) 



where 
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7(C) = (9 C + \/3 v/27 C2 + 256 C^) . (12) 
The analytical expression for r and 9 hy ^ and x therefore, given by the function u: 

r{f^,x) =u/x, (13) 



Cl 

and 



9{^,x) ~ arcsin-yu, 
where arcsin is defined as a continuous function of u with the range of [0, tt]. 
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III. MODIFIED DIPOLE COORDINATES 



The problem of the metric imbalance along field lines in the standard dipole coordinates originates from the power 
2 of the /z's denominator in eq. Therefore, one simple idea to reduce the steep metric distribution in the 
standard dipole coordinates shown in Fig.^a) is to use a coordinates (f) in which field-aligned coordinate /it', 

instead of /x, is defined as 



/.' = -^^^, for e<TTl2. (15) 
r 

It is easy to confirm that (/i',x, <j)) is also an orthogonal system. The metric of /i'-coordinate is given by 



/i^, = 1 / I V/i'l = 2 Vcos e /Q. (16) 

Hysell et al. P| used essentially the same coordinates as (/^', x, cf) for a plasma clouds simulation in midlatitude. (They 
used M = r/Vcos9 — — l//i', instead of /i'.) When one uses this coordinate system (/i', x, <t>), a care should be taken 
for the fact that the metric hf^i vanishes in the equator; h^i\^,^Q = 0. This was not a problem in the simulation by 
Hysell et al. ,4] since it was sufficient for them to use only a small part of the ^' space (—0.79 < fi' < —0.74). 

A simple remedy to avoid the singularity of ft.^' in the equator is not to place the grid point just on the equator. For 
example, the northern hemisphere is fully covered in practice by /i' < — e^', or 9 < tt/2 ~ eg, with a small positive 
buffer e^' or eg. Fig. O shows the grid points with iV^,/ x = 51 x 10 for the practically full range of the northern 
hemisphere by setting eg — 0.01. 

The transformation formula for the modified dipole coordinates (p' , x, 0) into the spherical coordinates are obtained 
by the same equation (jHJ by letting ( — {^x' /xY- 

Another form of modified orthogonal dipole coordinates proposed in this note is {if), Xi 4')^ where the new coordinate 
^p is defined through fi as 

ip = sinh"""^ (a fi) / a, (17) 

or its inverse transformation 

= sinh (a ■0) / a, (18) 

where a is a parameter that controls the grid distribution along dipole field lines, and a is defined as a = sinh^^ a. 
Note the identity of sinh^^ x — log (a; + Vl + x"^) ■ 
The metric of ■0 is given by 



= ^ ar^ cosh(aip) / (a&). (19) 

dtp 



The distribution as a function of V' when the control parameter a = 100 is shown in Fig.^b), which should be 
compared with ft,^ distribution shown in Fig.^a). It should be noted that the vertical scales in Fig-Qfa) and (b) are 
different for one order of magnitude. The basic idea that has lead to the transformation 117|) is to relax the steep 
gradient of the metric along /i in Fig.QJa) by local scale transformations. We want to make the fi grid spacing along 
the field line being "shrunk" only near the equator; /i ^ or ~ 0. Therefore, the grid "shrink rate" along the 
field line should be a function with a steeple-like peak at V' = 0- An example of such a function is 1/ cosh-0. The 
grid shrink rate is given by dip/dn since it is the reciprocal of the grid "stretch rate" dfj,/dip; see eq. (|19ll . Solving 
dip/diJ, = 1/ cosh-0, we get ^ = sinhip. 

In the limit of a 0, ip = ii. As the parameter a increases, grid points near the Earth (r — 1) along field lines, 
which are highly concentrated in the standard dipole coordinates (see the upper panel of Fig.|2J), move toward the 
equator along the field lines. The denominator d in eq. I|17|l is introduced to keep the transformed coordinate tp 
being always in the range of [—1, 1]. Fig. 0] shows grid points distribution when a — 100 for x = 101 x 10. 
The metric for the field line starting from 70°N on the equator is only 12 times larger than that on r = 1. The 
coordinate transformation by the sinh-function — applied to the cartesian coordinates — was also used in our numerical 
simulations of the magnetosphere ^ ,1Q.] . 

The coordinates transformations from {tp,x,4') into {r,9,(t)) are given by eqs. H13|l and (|14|l with eq. (|18|l . 

The relation between components of a vector v in the spherical coordinates {vr, vg, v^} and in the modified dipole 
coordinates {v^,v-^,V(j,} is given by the same form as that in the standard dipole coordinates: 

2 cos 9 / Q sin 
-sin6l/e 2coi 


The inverse transformation is given by the transverse matrix. 
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IV. DISCUSSION AND SUMMARY 



The standard dipole coordinates (/x, x, (j)) defined by eq. 1^ is not a good choice for a base grid in numerical studies 
since the metric contrast along each field line is too intense. Instead of the standard dipole coordinates, we propose 
to use the modified orthogonal dipole coordinates defined by 

, / sinh~^ (a cos srn^ 9 \ , , 

V smh a f J 

with a tuning parameter a of the metric distribution along field lines. The standard dipole coordinates is a special 
case in the limit of a 0. Fig.^] shows the case when a = 100 for the field lines with foot points located between 
65°N and 70°N. Since ■0's metric, or the grid distribution, is not sensitive to the change of a around this value, one 
would not have to perform its fine control. For other x range, pertinent a value will be easily found by visual checks 
of grid distribution images like Fig.^ 

For problems in which a symmetry around the equator between the northern and southern hemispheres is present, 
one can also try another form of modified orthogonal dipole coordinates defined by 



/ / ;N / \I^Ut>U sill U , \ 

{^i',xA)=\-- , ,0 , (22) 




in which /i' < — e, with a small positive buffer e for the northern hemisphere. 

Recently, a non-orthogonal dipole coordinates that is designed so that the lower-most constant-/^ surface coincides 
with a constant-r surface (i.e., a sphere) is presented 0]. It is straightforward and effective to apply the coordinate 
transformations presented in this note for that nonorthogonal dipole coordinates, too. 
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FIG. 1: The metric distribution along field the line starting from 70°N at r = 1. (a) /i^ for the standard dipole coordinates 
(b) for the modified dipole coordinates. Note that the vertical scales between the two panels (a) and (b) are different for 
one order of magnitude. 
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FIG. 2: The standard dipole coordinates {n, x, <P) = (— cos 9 / r^, sin^ / r, (f) in a meridian plane, = const. The grid points are 
distributed with equal spacings in each direction in the computational space and x- The total grid size is x iV^ = 101 x 10. 
There is no skip of grid points in the figure. The starting points of the field lines are between 65° N and 70° N at r = 1. The 
upper panel is a closer view. 
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FIG. 3: Modified orthogonal dipolc coordinates (/i', Xi 4') i^i ^ meridian plane. The coordinate fi' is defined as jj,' = — Vcos 9 / r, 
with /i' < — e^/, OT 9 < tt/2 — eg. The small buffer tfj^i or eg is introduced to avoid the vanishing metric in the equator, 
(/i^/ = at /i' = 0.) Here ee = 0.01. The total grid size is A'^/ x A^^^ = 51 x 10 in this "almost" northern hemispheric region. 



FIG. 4: Modified orthogonal dipole coordinates (tf), Xi 4') in meridian plane. The coordinate t/i is defined as tp = 
sinh~^ {a fi) / sinh^^ a. Total grid size is A'^^ x A^^ = 101 x 10. The control parameter a — 100 in this figure. Compare 
with the grid distribution of the standard dipole coordinates (Fig.|5J with the same grid size. 



